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We investigate the problem of an ultracold atomic gas in the superfluid phase flowing in the 
presence of a potential barrier or a periodic potential. We use a hydrodynamic scheme in the local 
density approximation (LDA) to obtain an analytic expression for the critical current as a function 
of the barrier height or the lattice intensity, which applies to both Bose and Fermi superfluids. In 
this scheme, the stationary flow becomes energetically unstable when the local superfluid velocity is 
equal to the local sound velocity at the point where the external potential is maximum. We compare 
this prediction with the results of the numerical solutions of the Gross-Pitaevskii and Bogoliubov-de 
Gennes equations. We discuss the role of long wavelength excitations in determining the critical 
velocity. Our results allow one to identify the different regimes of superfluid flow, namely, the LDA 
hydrodynamic regime, the regime of quantum effects beyond LDA for weak barriers and the regime 
of tunneling between weakly coupled superfluids for strong barriers. We finally discuss the relevance 
of these results in the context of current experiments with ultracold gases. 

PACS numbers: 03.75.Lm, 03.75.Ss 



I. INTRODUCTION 

The critical velocity of superfluid flow is one of the 
most fundamental issues in the physics of quantum flu- 
ids. A mechanism for the onset of dissipation is pro- 
vided by Landau instability [J 0, H, IH : if the excitation 
spectrum satisfies suitable criteria, there exists a critical 
flow velocity above which the kinetic energy of the su- 
perfluid can be dissipated via the creation of excitations. 
For uniform weakly interacting Bose gases, this critical 
velocity v c coincides with the velocity of sound waves. 
Several experiments have already been devoted to ob- 
serve and characterize the onset of dissipation in super- 
fluids made of ultracold atomic gases, including studies 
of Bose-Einstein condensates (BECs) in relative motion 
with respect to external perturbations, such as point-like 
defects [5[ or weak optical lattices @ . 

The present work was stimulated by the recent ex- 
periment by Miller et al. 0, where a one-dimensional 
(ID) optical lattice is used to measure the critical ve- 
locity of the superfluid flow of ultracold fermionic atoms 
in the BCS-BEC crossover. The lattice is produced by 
two intersecting laser beams focused in the central part 
of the trapped gas. A frequency difference between the 
two beams is used to produce a relative motion of the 
lattice and the gas. The number of fermion pairs which 
remain in the superfluid after the application of the mov- 
ing lattice is measured as a function of the relative ve- 
locity in order to determine the critical velocity for the 
onset of dissipation. The experiment confirms a theoret- 
ical prediction that superfluidity is most robust at uni- 
tarity [j| E3|- It also shows that the critical veloc- 
ity is very sensitive to the intensity of the lattice. This 
shares some similarity with the flow of a superfluid in 
the presence of a single potential barrier. In the case of 



ultracold fermions in the BCS-BEC crossover, the latter 
problem has been recently addressed in the framework of 
the Bogoliubov-de Gennes (BdG) theory by Spuntarelli 
et al. 11], who also found a strong dependence of the 
critical velocity on the barrier height. These works (see 
also, e.g., Refs. [H, EH EH EE EH) are recent fermionic 
examples of a wider field of investigations of critical cur- 
rent phenomena mainly studied, in the past, for Bose- 
Einstein condensed gases. We notice however that the 
similarities and differences between the behavior of the 
critical flow in the presence of a single barrier and of a 
periodic potential, as well as the role of quantum statis- 
tics have never been addressed in a clear and systematic 
way. 

Our work is aimed to establish an appropriate frame- 
work where one can compare the different situations 
(bosons vs. fermions and single barrier vs. optical lat- 
tice) and extract useful indications for available and/or 
feasible experiments. To this purpose we first introduce 
the hydrodynamic formalism in the local density approx- 
imation (LDA) and we assume a polytropic equation of 
state, which applies both to a Fermi superfluid at uni- 
tarity (i.e., when the s-wave scattering length a s is much 
larger than the interparticle distance, which in turn is 
much larger than the range of the interatomic potential) 
and to a Bose-Einstein condensate. In this way we derive 
a universal relation between the critical current for ener- 
getic instability and the maximum value of the external 
potential. We then compare the LDA hydrodynamic pre- 
diction with the results of the numerical solutions of the 
Gross-Pitaevskii (GP) and Bogoliubov-de Gennes equa- 
tions for bosons and fermions, respectively. This compar- 
ison allows us to identify different regimes of superfluid 
flow depending on the relevant energy and length scales 
of the system: a regime of hydrodynamic flow in LDA, a 
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regime of macroscopic flow beyond LDA, and a regime of 
weakly coupled superfluids separated by thin and strong 
barriers. From our analysis one can see that the BdG re- 
sults of Ref. for fermions at unitarity nicely fall into 
the LDA hydrodynamic regime. Conversely, the param- 
eters of the experiment in Ref. 0] are such that the LDA 
is not applicable, the healing length being of the same 
order of the lattice spacing. Indeed we find that the ex- 
perimental results are not reproduced by the LDA, as ex- 
pected, but we find that they also disagree with the BdG 
theory in a periodic potential with parameters similar to 
the ones of the experiment. This disagreement suggests 
that the interpretation of the experimental observations 
in terms of energetic instability of stationary superfluid 
flow in a uniform ID periodic potential remains an open 
issue. 



II. HYDRODYNAMICS IN THE LOCAL 
DENSITY APPROXIMATION 

Our starting point is a hydrodynamic theory in the lo- 
cal density approximation at zero temperature. Let us 
consider a superfluid which extends to infinity in three 
dimensions and is subject to a ID potential, V ex t(z), 
which has a finite maximum value V max . Let the su- 
perfluid be in a stationary configuration characterized 
by the density profile n(z) and the (quasi-) momentum 
P(z) along the z-direction. The LDA assumes that, 
locally, the system behaves like a uniform gas; thus 
the energy density e(n,P) can be written in the form 
e(n,P) = nP 2 /2m + e(n,0) and one can define the lo- 
cal chemical potential nip) = de(n, 0)/dn. The density 
profile of the gas at rest in the presence of the external 
potential can be obtained from the Thomas-Fermi rela- 
tion = ji{n(z)) + V cx t{z)- If the gas is flowing with 
a constant current density j — n(z)v(z), the Bernoulli 
equation for the stationary velocity field v(z) is 



* =2 



n(z) 



+ n{n) + V e xt(z), 



(1) 



where fij is the z-independent value of the chemical 
potential. This equation gives the density profile, for 
any given current j, once the equation of state /J-(n) 
of the uniform gas of density n is known. We assume 
the uniform gas to obey a polytropic equation of state, 



H{n) 



1 , and we consider two cases: i) a dilute Bosc 



gas with repulsive interaction, for which one has 7 = 1 
and a = g = AiTti?a s /m, where g is the interaction pa- 
rameter, m is the atom mass, and a s > is the s-wave 
scattering length; ii) a dilute Fermi gas at unitarity, for 
which u(n) = (1 + f3)E-p, where is a universal param- 
eter [13, Hll while E F = h 2 kl/(2m) and fc F = (S^n) 1 ^ 
are the Fermi energy and momentum, respectively, of a 
uniform non-interacting Fermi gas of density n, so that 
7 = 2/3 and a = (1 + (3)(3n 2 ) 2 / 3 h 2 /(2m). Using the 



equation of state one can write 



d 



mc s {z) = n—n(n) = jfi(n) . 



(2) 



where c s (z) is the local sound velocity, which depends 
on z through the density profile n(z). In a uniform 

gas of density n , the sound velocity is given by cf" 1 — 
[7^(n )/m] 1/2 . 

In the framework of LDA, the system becomes en- 
ergetically unstable when the local superfluid velocity, 
v(z) at some point z is equal to the local sound ve- 
locity, c s (z). If the external potential has a maximum 
at z = zq [i.e., V cx t{zo) — Knax], then at the same 
point the density is minimum, c s (z) is minimum and 
v(z) is maximum. This means that the superfluid be- 
comes first unstable precisely at z = zq. The condition 
for the occurrence of the instability can be written as 
m[j c /n c (z )] 2 = 7,u[n c (zo)] = janj(zo), where n c (z) is 
the density profile calculated at the critical current [T{|. 
By inserting this condition into the Bernoulli equation 
(H|), after a straightforward calculation one obtains the 
following implicit relation for the critical current: 
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It is worth noticing that this equation contains only z- 
independent quantities. It is also independent of the 
shape of the external potential: the only relevant param- 
eter being its maximum value Vmax- Moreover, it can 
be applied to both bosons and fermions. Its version for 
bosons in sl owly varying potentials was already discussed 
in Refs. [13, El (see also Ref. H2). 

In Figs. [1] and [H we plot the critical velocity ob- 
tained from the hydrodynamic expression ([3]) for bosons 
and fermions, respectively (thick solid lines). The up- 
per plots refer to the case of a single ID rectangu- 
lar barrier of width L and height Knax, the superfluid 
having an asymptotic constant density no and velocity 
vq = j/no at large distances from the barrier. The bot- 
tom plots refer to the case of a ID periodic potential 
V cxt {z) = KiaxSin 2 (q B z) = sE R sin 2 (q B z), where q B is 
the Bragg wave vector, _Er = ti 2 q 2 i /2m is the recoil en- 
ergy, s is the lattice strength, and the superfluid has an 
average density n and a macroscopic velocity defined by 
vo = j/no. In all cases, the critical velocity v c = j c /n-o is 
normalized to the value of the sound velocity in the uni- 
form gas, ci \ and is plotted as a function of V ms , x / ^j=o 
(see Appendix |A"|) . The limit V max / Hj=o —* corresponds 
to the usual Landau criterion for a uniform superfluid 
flow in the presence of a small external perturbation, i.e., 
a critical velocity equal to the sound velocity of the gas. 
In this hydrodynamic scheme, as mentioned before, the 
critical velocity decreases when Vmax increases mainly be- 
cause the density has a local depletion and the velocity 
has a corresponding local maximum, so that the Landau 
instability occurs earlier. When V mayi = fij c , the density 
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FIG. 1: (Color online) Critical velocity for energetic insta- 
bility of a dilute bosonic superfluid subject to a ID external 
potential. The critical velocity is given in units of the sound 
velocity of a uniform gas, ci°' , and is plotted as a function of 
the maximum of the external potential in units of the chemical 
potential /ij=o of the superfluid at rest. Top panel: the case of 
a single potential barrier. Bottom panel: the case of a periodic 
potential. Thick solid lines: prediction of the hydrodynamic 
theory within the LDA, calculated from Eq. ©. Symbols: 
results obtained from the numerical solution of the CP equa- 
tion @ for various values of L/£. Filled circles: L/£ = 0.5 
(gno/E-R, = 0.1); open squares: L/£ = 1 {gna/En = 0.4); 
filled triangles: L/£ = 1.57 (gno/En = 1); open trian- 
gles: L/£ — 3 (gno/En = 3.65); filled squares: L/£ = 5 
(gn /E R = 10); open circles: L/£ = 10 (gn /E R = 40). The 
thinner black solid lines are the tight-binding prediction @ 
for L/£ = 1 and 1.57. Dashed lines are guides to the eye. 
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FIG. 2: (Color online) Same as in Fig.[T]but for a superfluid 
of dilute fermions at unitarity. Thick solid lines: prediction of 
the LDA hydrodynamic theory, from Eq. <(3j . Symbols in the 
top panel: BdG results of Ref. [llj with LkF = 4. Symbols 
in the bottom panel represent our BdG results for a periodic 
lattice. The s / values are shown as follows: filled circles 
for Lk F = 0.497 {E F /E R = 0.1) and s = 0.1, 0.5, 1, 2.5, 3, 
and 3.5; open squares for Lk F = 0.886 (E f /Er = 0.3178) 
and s — 0.1, 0.5, 1, 1.5, 2.5, 3, and 4; filled triangles for 
Lk F = 1.11 {Ef/Er = 0.5) and s = 0.1, 0.5, 1, 1.5, and 2.5; 
open triangles for Lfcp = 1.57 (Ef/Er — 1) and s = 0.1, 0.5, 
1, 2.5, and 4; filled squares for LkF = 1-92 (Ef/Er = 1.5 and 
s = 0.5 and 1; open circle for Lfcp = 2.48 (Ef/Er = 2.5) 
and s — 2.5. The result for Lfcp = 0.497 and s — 0.1 being 
larger than unity is due to the numerical uncertainty, which is 
estimated to be of the order of 3% or less for all BdG points 
in the bottom panel of this figure. The thinner black solid 
lines are the tight-binding prediction ((9)1 for LkF = 0.497 and 
0.886. Dashed lines are guides to the eye. 



exactly vanishes below the barrier and hence the critical 
velocity goes to zero. 

The condition for the applicability of the LDA expres- 
sion ([3]) is that the typical length scale of the external 
potential must be much larger than the healing length 
£ of the superfluid. For a ID potential barrier of width 
L, this implies L/£ 1 [23|. In an optical lattice, the 
typical scale is the lattice spacing d = tt/qb] hi order to 
compare the two cases we define the barrier width L as 
half the lattice spacing, i.e., L = d/2 = 7t/(2q , b). The 
healing length of a dilute Bose superfluid of density no is 



£ = hj ' {2mgno) 1 / 2 , while for a Fermi gas at unitarity one 
has £ ~ which is consistent with the BCS coher- 

ence length £bcs = S^F/A gap , where vv — hk-p/rn and 
A gap is the pairing gap, which is of order E-p at unitar- 
ity. Effects beyond LDA become important when £ is of 
the same order or larger than L; they cause a smoothing 
of both density and velocity distributions, as well as the 
emergence of solitonic excitations, which are expected to 
play an important role in determining the critical veloc- 
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ity. This is well known for bosons, where a supercritical 
current results in the emission of shock waves in classi- 
cal hydrodynamics and of solitons in the Gross-Pitaevskii 
theory [H [H, HI • 

In closing this section, it is worth mentioning that 
within the LDA scheme one can easily calculate also the 
current ( velocity )-phase relation of the superfluid flow- 
ing through a single potential barrier. The phase <fi of 
the order parameter is related to the local velocity by 
v(z) = (h/M)W<f>, where M = m for Bose superfiuids 
and M — 2m for Fermi superfiuids. For a given asymp- 
totic velocity vo = j/no at z = ±00, one can directly 
integrate the phase across the barrier in order to ob- 
tain the phase difference 8(f). By subtracting the con- 
tribution of Vq, the phase difference can be defined as 
5(f) = {M/h) dz [v(z) — v }. In the LDA, the local ve- 
locity v{z) is constant under the rectangular barrier and 
one thus obtains 8(f) — ML[v(zq) — vo]/h. Examples are 
shown in Fig. [3] for bosons with L/£ = 10 (top panel) 
and unitary fermions with LUf = 4 (bottom panel). The 
predictions of the hydrodynamic theory within LDA are 
shown as solid lines for three values of the barrier height. 
Each line is drawn up to the value of 8(f> c for which the 
velocity becomes maximum. These lines are called stable 
branches in the current ( velocity )-phase diagram; on the 
right of the maximum there are no solutions within LDA 
(see the end of Section IIIII for details) . The maximum 
velocity coincides with the critical velocity v c plotted in 
the top panels of Figs. [T]and[2j 



III. GROSS-PITAEVSKII THEORY FOR 
BOSONS 

In order to account for effects beyond the LDA hydro- 
dynamic theory we use the GP theory for dilute bosons 
[2a , [271 [2a and the BdG equations for fermions at uni- 
tarity [29|. Let us first consider the case of a dilute 3D 
Bose-Einstein condensate at zero temperature subject to 
either a ID rectangular potential barrier or a periodic 
potential. The system is well described by the ID GP 
equation [28| : 
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FIG. 3: (Color online) Velocity-phase diagrams in the LDA 
hydrodynamic regime for a bosonic superfluid (top panel) and 
for a fermionic superfluid at unitarity (bottom panel) through 
a single barrier. The thick solid lines show the prediction of 
the LDA hydrodynamic theory. The symbols in the top panel 
are the results obtained from the GP equation; those in the 
bottom panel are taken from the BdG results of Ref. [Tlj . The 
barrier heights for bosonic case are Vm&x/ /J.j=a =0.05 (open 
circles), 0.2 (filled squares), and 0.4 (filled circles). Those 
for fermionic case are Kn ax / Hj=o = 0.042 (open circles), 0.17 
(filled squares), and 0.68 (filled circles). The barrier width is 
L/£ = 10 for the bosonic case, and Lkp = 4 for the fermionic 
case; these parameters are such that the corresponding curves 
for v c /c^ in the top panels of Figs. [1] and [2] are qualitatively 
similar. The dashed lines are guides to the eye. 
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where ^f(z) = \J n(z) exp[i</>(z)] is the complex order pa- 
rameter of the condensate, (f>(z) is its phase, and fi is the 
chemical potential. The local superfluid velocity is given 
by v(z) = (h/m)d z (f)(z) and the sound velocity in a uni- 
form condensate of density n is ci°^ = yj gn^/m. We 
numerically solve the GP equation in order to find the 
condition for the energetic stability of the condensate. 

In the case of a single barrier, we calculate the critical 
current by looking for solutions of the time-independent 
GP equation ([4]) at a fixed current j and for a given 
asymptotic density n 30]. We also use the time- 
dependent version of the GP equation in order to check 



the dynamical behavior of the unstable solutions. In 
the case of the periodic potential, the critical current 
is instead found by using the Bloch wave formalism and 
solving the linear stability problem for the GP energy 
functional as in Refs. [U, [32|, [33|. For the same case, 
we also use an alternative method, namely, a hydrody- 
namic analysis of the excitations, which is valid when 
the excitations that trigger the instability are those with 
wavelength much larger than the lattice spacing. For a 
superfluid with average density no and quasi-momentum 
P, this approach predicts the dispersion relation 
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where fuv and q are the energy and the wavenumber of 
the excitations. The energetic instability occurs when 

Emm 
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dn^dP 
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dn^dP 2 



(6) 



In practice, we calculate the energy density e(no,P) by 
solving the GP equation for a given average density no 
and quasi-momentum P of the superfluid. The value 
of P for which Eq. is satisfied is the critical quasi- 
momentum, P c . Finally, the critical velocity can be cal- 
culated by using the relation v c = (l/n )(de/dP)p c (see 
Appendix B for details). We have checked that the crit- 
ical velocity obtained in this way agrees with that ob- 
tained by means of the complete linear stability analy- 
sis within 1% in the whole range of gn^/E^ and V max 
considered in the present work. This confirms that the 
energetic instability in the periodic potential is driven by 
long-wavelength excitations, since the excitation energy 
of the sound mode is the smallest in this limit, as it oc- 
curs in the case of a single barrier. Equation © can be 
used also for unitary fermions, as we will see in Section 
IV. 

In Fig. [Tl we show the critical velocity in the case of 
the single rectangular barrier (top panel) and the peri- 
odic potential (bottom panel) for various values of L/£. 
In both cases, one clearly sees that the results of the GP 
equation approach the LDA prediction for L/£ 3> 1, as 
expected. The way of approach is, however, different. In 
the case of the periodic potential, v c exhibits a plateau 
for L/£ < 1 and small V max ; the plateau is instead absent 
in the case of the single barrier. The latter case is sim- 
pler, because the chemical potential and the excitation 
spectrum are fixed by the asymptotic density hq only 
and are unaffected by the external potential. In the limit 
of weak and thin barriers, our GP results agree with the 

analytic expression v c /cs — 1 — (3/4)(LI4iax/C5 n o) 2 ^ 3 
already derived in Ref. [2l| ; conversely, for thick barriers 
> 1), the Vmax -> limit of the LDA expression © 

is v c /c { s ] ~ 1 - V372(K iax /gno) 1/2 [25]. All curves ob- 
tained from the GP equation for the single barrier have a 
curvature which lies in between these two limiting cases. 
Instead, the periodic potential extends over the whole 
system and, consequently, for a given average density 
no, a change of V ma-X affects both /i and the excitation 
spectrum in a way that depends also on the value of the 
healing length £. In particular, if £ is larger than the 
lattice spacing and V max is not too large, the energy as- 
sociated with quantum pressure, which is proportional 
to l/£ 2 , acts against local deformations of the order pa- 
rameter and the latter remains almost unaffected by the 
modulation of the external potential. This is the regime 
of the plateau in Fig. [T] In terms of Eq. ©, this regime 
occurs when the left-hand-side is ~ P/m and the right- 



hand-side is 



„(0) 



so that the critical quasi-momentum 
obeys the relation P c /m — ci°\ which is the usual Lan- 
dau criterion for a uniform superfluid in the presence of 



weak impurities (see Appendix [B] for more details) . This 
regimes ends when, by increasing Vmax, the sound veloc- 
ity c s in the lattice [i.e., the right-hand-side of Eq. © 
evaluated at P = 0] becomes comparable or less than the 
maximum of the left-hand-side. Since, for small V max , 
H is proportional to no and the maximum of the left- 
hand-side of Eq. ([6]) is of the order of Os/m, the above 
condition is obtained for ji ~ Er. When Vmax is further 
increased, the chemical potential (i becomes larger than 
£7r, the density is forced to oscillate and v c /c^ starts 
decreasing . The system eventually reaches a regime of 
weakly coupled superfluids separated by strong barriers. 
This regime will be discussed in Section IVl 

It is worth noticing that, for Vmax < 2gno [351 ] . loops 
("swallow tails") appear in the Bloch band structure of 
the superfluid in the periodic potential as a consequence 
of the nonlinearity of the GP equation. These states and 
their stability have been deeply discussed in Ref. [32j and 
are accounted for in our calculations. The sets of points 
in the bottom panel of Fig. [T] which are closest to the 
LDA curve (i.e., for L/£ = 10 and V max / fij = o < 0.9, or 
for L/£ = 5 and V max /lij=o < 0.8) fall into this regime: 
the critical velocity is determined by the energetic in- 
stability occurring along the swallow tails and the crit- 
ical quasi-momentum is outside the first Brillouin zone, 
P c > hqB- Our numerical results for P c well agree with 
those of Ref. [13]. 

Finally, from the solution of the GP equation one can 
calculate the current (velocity )-phase relation of the su- 
perfluid flowing through a single potential barrier. Ex- 
amples are given in the top panel of Fig. [3] The figure 
shows that, for L/£ = 10, the results of the GP calcu- 
lations are close to those in the LDA limit not only for 
the critical velocity (the maximum value of v) but also 
for the shape of the curve on the left of the maximum 
(stable branch). On the right of the maximum (unsta- 
ble branch) the GP theory still works, providing solitonic 
solutions which are instead not accessible to the hydro- 
dynamic theory without dispersion. We note that, when 
approaching the maximum from the right, the amplitude 
of solitons predicted by the GP equations vanishes and 
the two, stable and unstable, branches merge together 



IV. BOGOLIUBOV-DE GENNES THEORY FOR 
UNITARY FERMIONS 

A similar analysis can be done for a superfluid Fermi 
gas in the BCS-BEC crossover at zero temperature. A 
suitable approach consists in the numerical solution of 
the BdG equations for the fermionic quasiparticlc ampli- 
tudes u;, and Vi 12911 . 



H'(t) A(r) 
A*(r) -JT'(r) 



Ui(r) 
Vi(r) 



Ui(r) 
Vi(r) 



(7) 



with H'(r) = —h 2 V 2 /2m + V cxt — \x. These equations 
can be used to calculate the order parameter A(r) = 
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— g^2i u i{ r ) v i ( r )> the chemical potential fi, and the en- 
ergy density e. 

For a superfluid Fermi gas in the presence of a single 
rectangular barrier, the BdG equations have already been 
used by Spuntarclli et al. [1 11 ] to calculate the critical 
current in the case of Lkp = 4. Their results at unitarity 
are shown in the top panel of Fig. [2] (dots) where they 
are compared with the LDA prediction ([3]). The figure 
shows that, for this value of Lkp, the BdG results are 
indeed close to the LDA curve; in this regime, also for 
fcrmions, the critical velocity is mostly determined by 
the local depletion of the density where the potential is 
maximum. Concerning the LDA curve, we notice that 
in the Knax - * limit one can easily obtain the analytic 
expression 

v c /cf ~ 1 - V2 [V max /(1 + p)E F ] 1/2 . (8) 

In the lower panel of Fig. [2] we instead show our BdG 
results for unitary fermions in a periodic potential for var- 
ious values of Lfcp- These results are obtained by using 
Bloch functions and solving the BdG equations with the 
same self-consistent procedure of Ref. !36|. In particular, 
the energy density e(nn, P) is calculated as in Eq. (8) 
of Ref. [36| (see also [33] )■ From the energy density we 
then obtain the critical quasi-momentum and velocity by 
means of Eq. ([6]). 

A striking similarity emerges from Figs. [1] and [2) re- 
vealing that bosons and unitary fermions have a rather 
similar behavior. The rapid approach of the BdG results 
to the LDA curve for Lkp > 1 in the lattice is consistent 
with the case of a single barrier, as well as with the behav- 
ior of bosons in the corresponding regime, L/£ > 1. Our 
BdG calculations, which include an analysis of the quasi- 
particle energy spectra, show that the role of fermionic 
pair-breaking excitations, not included in Eq. ((6|), is neg- 
ligible at unitarity, except for small Vmax and for very 
low densities, such that Lfcp <C 1 [38 |. These fermionic 
excitations are instead expected to become important on 
the BCS side of the BCS-BEC crossover [M Cl- 
in the bottom panel of Fig. [3] we show the BdG results 
of Ref. [IH for the current (velocity)-phase relation of the 
unitary Fermi superfluid flowing through a single poten- 
tial barrier. As in the case of bosons (top panels) one 
can see that, for Lk? larger than 1, the BdG results are 
close to the LDA limit both for the shape of the stable 
branch and the height of the maximum. 



V. WEAKLY COUPLED SUPERFLUIDS 

By increasing the height of the barriers, such that 
V max ^S> n, the critical velocity decreases. In the LDA 
regime (L/£ ^> 1 or Lkp ^> 1) this decrease is fast be- 
cause the superfluid density locally follows the shape of 
the external potential and rapidly vanishes under a bar- 
rier as soon as V max exceeds [i. Conversely, when £ is 
of order L or larger, the coupling between the superfluid 



regions on the two sides of a barrier remains significant 
also for larger values of V max , leading to the physics of 
macroscopic tunneling and Josephson junctions in weakly 
coupled superfluids (for bosons, see for instance [39j]). 

In this regime, one finds the well-known characteristic 
sinusoidal current-phase relation, i.e., j(S<p) — j c sin(<50), 
where 8(f) is the phase difference of the superfluid order 
parameter across the barrier [40(. In the periodic lattice 
we recover the tight-binding expression for the energy 
density, which is given by e(n ,P) = e(n o ,0) + 8j[l — 
cos(7rP/P c dg C )]- Here, P d go is the quasi-momentum at 
the edge of the Brillouin zone, namely, P e dge = ^<7b for 
BECs of bosonic atoms and P e d go = %b/2 for superfluids 
of fermionic atoms. The quantity 8j — noP c 2 dge /7r 2 m* 
corresponds to the half-width of the lowest Bloch band 
and is related to the tunneling energy associated with the 
Josephson current. The effective mass m* is defined via 
the relation 1/m* = (l/n )d 2 P e(n , P)\ P=0 0,53- 

A consequence of the sinusoidal shape of the energy 
density in the tight-binding limit is that the critical quasi- 
momentum approaches the value P c = P dge/2 and the 
critical velocity takes the following simple dependence on 
the effective mass: 



When s — V max /E K ^> 1, the critical velocity decreases 
according to m/m* oc exp(— 2y/s). The values of v c ob- 
tained from Eq. ©, with m* extracted from the GP cal- 
culation of e(no,P), are plotted in the lower panel of 
Fig. Q] (thinner black solid lines) for L/£ = 1 and 1.57 in 
the region of V max / fJ>j=o > 2. The agreement with the 
results of the complete stability analysis is remarkable. A 
similar agreement is also obtained for unitary fermions, 
when e(no,P) is taken from the BdG calculations, as 
shown in the lower panel of Fig. [2] 

VI. FINAL REMARKS 

In all panels of Figs. Q] and [2] one can distinguish three 
limiting cases: i) a regime of hydrodynamic flow in the 
local density approximation for large L/£ (or Lkp), cor- 
responding to the points close to the thick solid lines; ii) 
a regime of macroscopic flow through thin and weak bar- 
riers, where the LDA is not applicable, i.e., for L/£ < 1 
(or Lk F < 1) and V ma ^/fi < 1; iii) a regime of weakly 
coupled superfluids separated by thin and strong barri- 
ers, i.e., for L/£ < 1 (or Lkp < 1) but V max /fJ> 1. 

Concerning bosons, the experiments performed in 
Ref. [iH with condensates in deep optical lattices fall 
into the third regime and were consistently interpreted 
in terms of the dynamics of a chain of Josephson junc- 
tions. An investigation of energetic and dynamical insta- 
bilities of bosons in shallow lattices was instead reported 
by the same group in Ref. 0], in a range of parameters 
corresponding to the second regime, namely where the 
critical velocity for the energetic instability is close to 



7 



the sound velocity of the uniform superfiuid, as in the 
standard Landau criterion. In this experiment, as in the 
one of Ref. p| , a sharp determination of the Landau criti- 
cal velocity was, however, hindered by the inhomogcncity 
of the trapped condensate. 
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FIG. 4: (Color online) Comparison between theoretical re- 
sults and the experimental data for the critical velocity of 
unitary fermions in an optical lattice. The filled triangles are 
our BdG results calculated for Ef/Eh — 0.5. The solid line 
is the prediction of hydrodynamics in LDA. The open squares 
are the experimental result of Ref. 0, normalized by using 



the density at the e 
to E f /Er = 0.56. 



waist of the lattice, which corresponds 



the BdG theory (e.g., the discrepancy between the mean- 
field result and the Monte-Carlo result of (3) is factorized 
out by the choice of the normalization in Fig. [H This 
suggests that other effects, including the inhomogeneity 
of both the density and the lattice intensity, as well as 
the nonstationarity of the process, can be important in 
explaining the experimental observations. These effects 
could be better understood by repeating the same type 
of experiment with bosons, where the GP theory is reli- 
able, full 3D GP simulations are feasible and a number 
of theoretical approaches are available to treat the case 
of time-dependent optical lattices (as, for instance, in 
Bragg spectroscopy experiments), from linear response 
theory to full nonlinear dynamics. 

In view of the importance that energetic instabilities 
have in the description of superfiuid phenomena, fur- 
ther experimental and theoretical investigations on the 
behavior of ultracold atomic superfluids in optical lat- 
tices and in the presence of obstacles are worth pursuing. 
An important step is the choice of a suitable geometry. 
Elongated systems with strong transverse confinement 
are good candidates, because their excitation spectrum 
is simpler. A new interesting perspective is also pro- 
vided by the availability of toroidal confining potentials 
for ultracold gases 0, |4^], in which one can produce a 
stationary superflow and explore the dissipation induced 
by an external perturbation along the torus [46[. 



The critical velocity of fermions in an optical lattice 
is measured in Ref. The effects of the inhomogene- 
ity of the density is reduced by producing the lattice in 
the central part of the atomic cloud. However, the lat- 
tice itself turns out to be inhomogeneous. This causes 
some uncertainties in the determination of both kp and 
VJnax, which are needed for the comparison with a the- 
ory for a uniform superfiuid in a uniform lattice. The 
waist of the two laser beams is smaller, but not much 
smaller than the Thomas-Fermi radii of the trapped gas. 
A reasonable choice consists of taking the density no as 
the density measured at the e -2 waist of the lattice, as 
suggested in Ref. 0. This yields E F /E R ~ 0.56, which 
means Lkp ~ 1.1. In Fig. Uwe show the experimen- 
tal data (open squares) for this choice of kp. Our BdG 
results for Ep/E R — 0.5 are shown as filled triangles, 
while the LDA result is represented by the thick solid 
line. Here v c and Kn ax are normalized to the sound ve- 
locity ci -* = + (3)/2>vp and the chemical potential 
(1 + 0)Ep of the uniform gas, respectively [18| . The be- 
havior of the experimental data significantly differs from 
the theoretical predictions. The disagreement remains 
qualitatively the same even for different choices of kp and 
V max within the experimental uncertainties [43[ . The fail- 
ure of the LDA is expected, being consistent with the fact 
that here Lkp is of the order of unity. The disagreement 
with the BdG calculations is instead puzzling, especially 
if one considers that most of the expected inaccuracy of 
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APPENDIX A: CRITICAL VELOCITY IN THE 
LDA HYDRODYNAMIC THEORY 

In order to plot the LDA results in Figs. Q] and O we 
first rewrite Eq. |T]) in terms of the velocity vo = j/no- 
Its critical value, v c , is determined by the condition 
that the local superfiuid velocity v(z) coincides with 
the local sound velocity c s (z) [Eq. at the position 
z = Zq, where the potential takes a maximum value 
K xt(zo) = Vma x- Using this condition, v(z ) = c s (z ) = 
y/a , yn' y (zo)/m, together with the continuity equation, 
n(zo)v(zo) = novo, we write Eq. fl} as 

^=(i+^)(ir + V'<->+^, (ad 
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(mv 2 /2)/ariQ, and V m 



where /2 C = fj,j c /wig, v t 
V ma x/an^. 

In the case of a single potential barrier, one finds jl c = 
v 2 + 1 , so that Eq. (|A1|) becomes a third order algebraic 
equation for v 2 ^ 3 when 7=1 (bosons), and a fourth order 

1/2 

equation for v c when 7 = 2/3 (unitary fermions). Once 
this algebraic equation is solved, one can use (J,j=o = anj 
[i.e., /ij=o = gno for bosons and fJ,j=o = (1 + 0)E-p for 
unitary fermions] and plot v c as a function of Vmax/M^o- 
In the case of the periodic potential the procedure is 
similar, but the chemical potential must be calculated nu- 
merically for each value of the average density tlq . Equiv- 
alently, if Z\ is a position where V CY x(z\) = 0, one has 
jl c = n^ 2 v 2 + nj, with h\ = n(zi)/riQ. We thus ob- 
tain an equation for v c similar to the one of the single 
barrier, but where h\ has to be determined numerically. 
Using fij—o obtained by a separate calculation, we can 
parametrically plot v c as a function of V max / jJ,j=o- 



APPENDIX B: DETERMINATION OF THE 
CRITICAL QUASI-MOMENTUM IN A LATTICE 

In this Appendix, we explain how we calculate the crit- 
ical quasi- momentum P c from the relation ([6]) . As an ex- 
ample, in Fig. [5l we show the case of the unitary Fermi 
gas at Ef / En = 1, where P c is in the first Brillouin zone 
[471 ]. Here, we drop the subscript "0" of the average den- 
sity no for simplicity. 

For the uniform case (s = 0), the energy density 
can be written as e(n, P) = nP 2 /2m + e(n, 0), and 
thus d 2 e(n,P) — d 2 e(n 7 0) = d n ^{n) — l/(nn) and 
dpe(n,P) = n/m. This implies 



d 2 e d 2 e 
dn 2 ~dP 2 ~ 
2 e 



_ r (o) 



dndp 



P 

m 



(Bl) 
(B2) 

where k is the compressibility and cf^ is the sound 
velocity of the uniform system: ci°^ = ^(1 + /3)/3vf 
(~ 0A4vp from the mean-field theory, and ~ 0.37«f 
from Monte Carlo calculations) for unitary fermions and 

ci -* = yj gn/m for bosons. These results are shown as 
blue and red dashed lines in Fig. [5] Their crossing point, 
P/m = Cs , gives P c . This procedure is consistent with 
the Landau criterion in a uniform gas. 

For the superfluid in a lattice (s ^ 0), since dpe = 
at the edge of the Brillouin zone, d n dpe is zero there and 
has a maximum in the first Brillouin zone (see the red 
solid line in Fig. [5|). For the same reason, there exists 
a new point P c ,dyn at which d p e — and the curve for 
y/d^edpe is bent downward (see the blue solid line in 
Fig.©. ThisP c> dyn, which does not exist for the uniform 
case, is the critical quasi-momentum for the dynamical 
instability of long-wavelength excitations; the excitation 
energy fku becomes complex for P > P c , dyn- The critical 




P/P 



edge 



FIG. 5: (Color online) Schematic picture for the determi- 
nation of the critical quasi-momentum P c for the energetic 
instability using the hydrodynamic relation (|6]). Here we take 
the unitary Fermi gas at Ef/E-r, — 1 as an example. Dashed 
lines correspond to the uniform case (s = 0) while solid lines 
to s — 1. The value of P c is given by the abscissa of the cross- 
ing point between the curve representing d n dpe (red line) and 
(dledpe) 1/2 (blue line). The critical quasi-momentum P c , dyn 
for the dynamical instability of long-wavelength excitations is 
given by a zero of (<9^e d%e) 1 ^ 2 , where dpe = 0. 



quasi-momentum for the energetic instability P c is still 
given by the crossing point of d n dpe and y/d 2 edpe (red 
and blue lines) and thus P c < P Cj dyn- 



Note that the value of \Jd\eSpe (blue curve) at P = 
coincides with the sound velocity c s in the lattice, which 
is decreasing function of s (48J. On the other hand the 
slope of d n dpe (red curve) is of order \/m* , where m* is 
the effective mass. Since c s cx 1/y/m* the slope of d n dpe 
approaches zero faster than c s with increasing s. For 
large s, the red curve takes a sinusoidal shape and the 
crossing between the two lines occurs closer and closer to 
P = P c dg /2. In the same limit one finds P c = P C; dyn- 

It is also worth noticing that, for both bosons and 
fermions in the lattice, there exists a value of average 
density n such that the critical momentum in the uni- 
form case, mcs , is equal to the critical quasi-momentum 
in the tight-binding limit, P d g o/2. In this case, we 
have numerically checked that P c is almost independent 
of the lattice strength s in the whole range s > 0. 
With fermions, this condition occurs when Ef/Er = 
(3/16)(l + /3) _1 ; for bosons, it occurs when gn /pR = 
l/y/2. For smaller (larger) density, the critical quasi- 
momentum P c is a monotonic function of s approaching 
the asymptotic value P e d ge /2 from below (above). 
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and we use the following fitting function: 
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FIG. 6: (Color online) Energy density for fermions at uni- 
tarity as a function of the quasi-momentum P for various 
values of the lattice strength s. Here we take Ef/E-r = 1 
(kpL = 1.57) as an example. The inset shows the magnifica- 
tion of the result for s = 2.5; the dotted line is a sinusoidal 
curve with the same amplitude. 



APPENDIX C: BAND STRUCTURE OF THE 
SUPERFLUID UNITARY FERMI GASES IN 
PERIODIC POTENTIALS 

For superfluid Fermi gases in an optical lattice, solving 
the BdG equations with a large energy cutoff Eq is rather 
time-consuming. In the present work, in order to obtain 
the P-dependence of e and fi for various parameters, we 
calculate them for a limited number of points (typically 
seven or more) in the first Brillouin zone for each Ef / Pr 



f(n ,P) =/(n 0) 0) 



+ 5(no)i - I 1-cos 



P 



P. 



edge 



with the band width, 



S(n ) = f(no, Podge) - f(n o ,0). 



1/7, 

) 

(CI) 
(C2) 



The shape parameter rj can take all values form 1 to 00: 
the case r\ = 1 corresponds to the sinusoidal form for 
strong lattices, while 77 = 00 reproduces the quadratic 
energy dispersion of the translationally invariant system 
without lattice (s = 0). Note that dpf(n ,P) = at 
P = and P e dge like e and \i (here, we do not consider 
situations where the band structure exhibits "swallow 
tails"). We find that the above fitting function (|C1[) very 
well reproduces the numerical BdG results of e and fj, ex- 
cept for the cases of very small s. Even for s = 0.1, which 
is the smallest non-zero value studied in the present work, 
the relative error is typically less than 0.1%. 

In Fig.[6l we show the energy density e(n, P) of the uni- 
tary Fermi superfluid in a periodic potential calculated by 
the BdG equations for E-p /Er = 1 and for three values of 
the lattice strength, s = 0.1, 1, and 2.5. The correspond- 
ing values of the shape parameter rj are i] = 5.05, 1.68, 
and 1.13. Note that the P-dependence of e is rather close 
to a quadratic form for the weakest lattice (s = 0.1) and 
is almost sinusoidal for the strongest lattice (s = 2.5), as 
one can see in the inset of Fig. [5] 
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